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Abstract 


Sixth-order compact finite difference method is presented for solving the 
one-dimensional KdV-Burger equation. First, the given solution domain 
is discretized using a uniform discretization grid point in a spatial direc- 
tion. Then, using the Taylor series expansion, we obtain a higher-order 
finite difference discretization of the KdV-Burger equation involving spa- 
tial variables and produce a system of nonlinear ordinary differential equa- 
tions. Then, the obtained system of a differential equation is solved by 
using the fourth-order Runge-Kutta method. To validate the applicability 
of proposed techniques, four model examples are considered. The stabil- 
ity and convergent analysis of the present method is worked by using von 
Neumann stability analysis techniques by supporting the theoretical and 
mathematical statements in order to verify the accuracy of the present so- 
lution. The quality of the attending method has been shown in the sense 
of root mean square error Lz and point-wise maximum absolute error Loo. 
This is used to show, how the present method approximates the exact so- 
lution very well and how it is quite efficient and practically well suited for 
solving the KdV-Burger equation. Numerical results of considered exam- 
ples are presented in terms of Lz and Log in tables. The graph of obtained 
present numerical and its exact solution are also presented in this paper. 
The present approximate numeric solvent in the table and graph shows 
that the numerical solutions are in good agreement with the exact solution 
of the given model problem. Hence the technique is reliable and capable 
for solving the one-dimensional KdV-Burger equation. 
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1 Introduction 


Partial differential equations (PDEs) are the mathematical equations that 
are significant in modeling physical phenomena that occur in nature. Ap- 
plications of PDEs can be found in physics, engineering, mathematics, and 
finance. Examples include modeling mechanical vibration, heat, sound vi- 
bration, elasticity, and fluid dynamics [4]. Some applications in real life 
phenomena appear as nonlinear fractional order of PDEs. The fractional 
nonlinear PDE describes a variety of important physical phenomena, acous- 
tics, elector-chemistry, materials science, fluid dynamics[2], optics, and visco- 
elasticity [1]. The other type of nonlinear PDEs is also Korteweg-de Vries 
Burgers’ equation. Korteweg-de Vries Burgers’ equation is a model mathe- 
matical statement for a wide class of nonlinear wave models of fluid in an 
elastic tube, liquid with small bubbles and turbulence [21]. This Burgers- 
Korteweg-de Vries (Burgers-KdV) equation has wide applications in physics, 
engineering, and fluid mechanics. The Poincaré phase plane analysis reveals 
that the Burgers-KdV equation has neither nontrivial bell-profile travelling 
solitary waves nor periodic waves [9]. 

The Burgers-KdV equation is the simplest form of the wave equation 
containing the nonlinearity, in which dispersion and the dissipation term 
occur [9]. It arises from many physical contexts, for example, the propagation 
of glandular bores in shallow water [5, 18], the flow of liquids containing 
gas bubbles [13], the propagation of waves in an elastic tube filled with a 
viscous fluid [36], weakly nonlinear plasma waves with certain dissipate effects 
[13, 16], and the cascading down the process of turbulence. It is also widely 
used as a nonlinear governing model in the crystal lattice theory, the nonlinear 
circuit hypothesis, and atmospheric dynamics [10]. 

Travelling waves or solitons as solutions to the Korteweg-de Vries equation 
(KdV), which is a nonlinear partial differential of mathematical statement 
four-third-order, have been of interest already for 150 years [37]. The Kdv 
equation was discovered in 1895 by Korteweg and de Vries [23], but this 
equation was forgotten for a long time. So that, the KdV equation is widely 
recognized as a paradigm for the description of weakly nonlinear long waves 
in many branches of physics and engineering. It describes how waves evolve 
under the competing but comparable effects of weak nonlinearity and dy- 
namic dispersion [37]. Random waves are an important subject of haphazard 
PDEs [11]. The nonlinear problems are solved easily and elegantly without 
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converting to linear form the problem by using the Adomian decomposition 
method [22]. A feature common to all the traditional methods is that they 
are using the transformations to reduce the equation into a more simple 
mathematical statement and then solve it. Unlike classical techniques, the 
nonlinear equations are solved easily and elegantly without transforming the 
equation by using the Adomian decomposition method [21]. This wonder- 
ful technique has many advantages over classical techniques [22]. Mainly, it 
avoids the process of converting linear form and perturbation to find solutions 
to a given nonlinear equation. It provides an efficient explicit solution with 
high accuracy, minimal calculation, and avoidance of physically unrealistic 
assumptions [21]. 

Many authors studied the exact and numerical solutions of the KdV equa- 
tion. Wadati [34] first introduced and studied the stochastic KdV equations 
and gave the diffusion of solutions to the KdV equation under Gaussian noise. 
Xie firstly studied research on Wick-type stochastic KdV equation on white 
noise spaces and showed the auto-Backlund transformation and the exact 
white noise functional solutions in [16]. Furthermore, Chen and Xie [6, 8] and 
Xie [7, 38] worked on some Wick-type stochastic wave equations using the 
white noise analysis method. Recently, Ugurlu and Kaya [14] gave the tanh 
function method, Wazzan [35] showed the modified tanh-coth techniques, 
and these methods have been applied to derive nonlinear transformations 
and exact solutions of nonlinear PDEs in mathematical physics [18]. In those 
driving analytical methods, there is time consumed and very tired because it 
requires long computational time. To save our time and obtain an accurate 
approximate solution of the KdV equation, we use numerical methods. 

Different numeral techniques give accurate denotative solutions for the 
KdV equation. To obtain the numerical solution of the KdV equations, the 
same researcher was used the traditionally, mesh-dependent methods such 
as the finite-difference method, finite element techniques, and boundary el- 
ements method. However, these methods have a slow rate of convergence, 
instability, low accuracy, and difficulty of implementation in complex ge- 
ometries [4]. Even though the accuracy of the aforementioned methods is 
promising, they require large memory and long computational time. Besides 
funding, the methods are not suitable for higher-dimensional and problems 
involving complex geometries. Hence, the treatment of mesh-size presents 
several difficulties that have to be addressed to ensure the accuracy of the 
solution, and efficiency of the method applied [4]. Indeed scientists in the 
field of computational mathematics have been trying to develop numerical 
techniques by using computers for further application of different model prob- 
lems. From those techniques, higher-order difference schemes and differential 
quadrature methods are more powerful techniques. 

In many application areas, such as aeroacoustics [26] and electromagnetic 
[29], the propagation of acoustic and magnetic force waves needs to be accu- 
rately simulated over very long periods of time and far distances. Especially, 
the nonlinear partial differential equation is derived that plays a role in a 
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normal form, that is, in the first approximation, it determines the behavior 
of all solutions for the original boundary value problem with initial condi- 
tions from a sufficiently small neighborhood of equilibrium [20]. Therefore, in 
order to reduce the accumulation of errors during this process, the numerical 
algorithm must be highly accurate. To accomplish this goal, high-order com- 
pact finite difference schemes have been developed for simulating the solution 
of applications to wave propagation like [31, 32] and surface water modeling 
(27, 28]. The investigation of numerical solutions for nonlinear partial differ- 
ential equations plays an important role in mathematics, physics, and other 
applied science areas. 


High-order finite difference schemes can be classified into two main cate- 
gories: explicit schemes and Pade-type or compact schemes. These hardcore 
techniques compute the numerical derivatives directly at each grid by us- 
ing large stencils, while compact schemes obtain all the numerical derivatives 
along a grid line using smaller stencils for solving linear and nonlinear systems 
of equations. Experience has shown that compact schemes are much more 
accurate than the corresponding explicit scheme of the same order. Hence 
motivating this performance of implicit (Pade-type) scheme, this paper aims 
to construct an efficient and accurate sixth-order compact finite difference 
numerical method for solving a one-dimensional KdV-Burger equation in the 
application of wave propagation. 


Statement of the problem 


Consider a homogeneous KdV-Burger equations considered in [24] given by 
U,+auvu, + weer = 9, (a, t) € (0,b) x (0,T), (1) 

with initial and boundary conditions, respectively, 
U(a,0) = Uo(x), 0<a<b, (2) 


U(0,t) = U,(t),U(b,t)=U,(t), O<a<T, (3) 


where U,(t) are U,(t) are periodic smooth functions and 0 < T < oo. Now 
we define a mesh size h and At as constant grid points. They are defined 
by drawing distant horizontal and vertical lines of distance “h” and “At”, 
respectively, in “x” and “t” directions. These lines are called grid lines, and 
the point at which they interact is known as the mesh point. This fabric 
topology point that lies at end of the domains, is called the boundary point. 
Points that lie inside the region are called interior points. The goal is to 
find the approximate solution of “U(x, t)” at the interior mesh points and to 
investigate how it converges too. Thus to do this, first, we discretized the 
solution domain concerning spatial direction as 
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where @j;41 = 2; +h and j =1(1)M. Moreover, M is the maximum number 
of grid points in the z-direction. Thus, to accomplish our work, the present 
paper is organized as follows. Section 2 is a description of numerical methods. 
Sections 3 is stability, consistence, and convergence analysis. Section 4 is the 
results of numerical experiments. Sections 5 is the discussion, and Section 6 
is the Conclusion. 


2 Description of proposed numerical method 


2.1 Spatial discretization 


In the high-order compact finite difference methods, the spatial derivatives 
in the governing PDEs are not approximated directly by some finite differ- 
ences. They are evaluated by some compact difference schemes. We assume 
that U(a,t) is a sufficiently smooth function and has continuous higher-order 
partial derivative on its domain. For the sake of simplicity, let U(a;,t) = U; 
OPU(ax;,t 
and let “Us. 
be 

of U(a,t) concerning with spatial variable x. By using the Taylor series 
expansion, we have 


= O8U; in which p > 1, be pth-order partial derivative 


h? a» ho os ht a4 
Uj41 = Uj; + hou; + gy aU + yy Oe; + qi t+ ; (4) 
h? h3 h4 
Uj-1 =U; — hou; ia oy O2U; 3 O3U; Al OtU; aan (5) 
ARP oe Be a 28 yO a2 
4h? = 8h3 5 16h4 
Uj-2 =U; — 2hOU; + 5-020; — GZ -aU; + O20; (7) 
Now subtracting (7) from (11), we obtain 
OU; = Uj41 — Oj-1 h? OU; ht OU; SS (8) 
: 2h GO fs 


Hence from (10), we obtain the second-order central finite difference equation 
for the first-order partial derivative of U(a,t) concerning spatial variable x 
given by 

Uj+i ~ Uj-1 


bp Uj = oh 


+ Ti, (9) 
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h? 
where T; = ——03U; = | O2U (2;,#) is its local truncation error. Again 


subtracting (9) from (8), we obtain 


8h3 ., 32h 128h" 
U; = U;25 — U;_5 — 4hOU; 5 
ee eR anu 5040 


OFU; —-++. (10) 


Now substituting (10) into (10), we obtain 


Bh 98 gti Uj-1 _ h* o3,, |) _ 82h? a5, 
BU; = Ujy2—-Uj—2 n(4 a 5 aU; ) op oe 


Collecting like terms in the above difference equation, we obtain 


Ujse-— 2(Uj41 — Uj) —Ujz22 TH 31h4 
au, = + j j J OU; atu;—--». 
gad 2h3 60 7 7% 2520” 1) 


Hence from (14), we obtain the second-order central finite difference equation 
for third-order partial derivative of U(«,t) concerning spatial variable x, given 
by 

Uj42 — 2(Uj41 — Uj—1) — Uj~2 


3 — | 
O.U; — 2h3 t To, (12) 
(i eee . ; : ere 
where Tz = ———0?U; is their local truncation error. Again substituting 
(11) and (7) into (9), we obtain 
h? 2» he as 
U; + h0,U; + —02U; + —O02U;4+--- 
sOu, = 2! 3! 
cx J ah 
h? h3 
@ hou; 4 1 82U; 31 03U; -) - 
(2h a 1; 
h? ns 
6\DU; =0,U; aU; I" ou, “pty, ! O8U; +++ 
cmc ce Rca 2 - 0 iz 362880 ” ; 
h2 
6\9U; =0,U; + —O3U; + A osu, T: 13 
cx! J. g 6 xj 120 x + 3) ( ) 
hs 
where T3 = ——0°U; is their principal local traction error. Again substitut- 


ing (11)-(9) into (15), we obtain 


(U; + 2h0,U; +--+) —2((U; + hd,U; +--+) — (Uj; — hOU; —---)) 


Qy. = 
oe U; 7 2h3 


(U; PhO 4325) 


One + T» 
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8h 
Qy. ea 307, 4 5U, 4 77T. 9 
ce Onn mn v ORU 315 Oe a aot ce 


Thus from this, we obtain the second-order difference equation given by 


Qu; = BU; °U; 4 U; +T. 14 
ocx U; 0, Uj; ie a h asu U; ao U; 4 4) ( ) 
BLP pee 3 ae : 
where Ty = 7976 eV is its principal local truncation error. Now using (2) 


and applying the successive derivative on it in terms of spatial variable z, 

we obtain the fifth- and sixth-order central difference scheme, and they are 
given by 

PU; = 82 (500, (15) 

Cet (sO (oe). (16) 


Now substituting (15) into (13), we obtain 


Ges 3 (3) (5(2) 
su ae Pag ee La (08 U;) +Ts, 
h hd 
Oy, ay, 4 Ou. 4 (3) (sq 
5DU; = aU; + IOU; + 58 (08 Uj) +Ts, 
OU5 gy, gay. Pgs (5(2) 
Ae = QU; — 8 U; — —5@) (62U;) - Ts. (17) 


Again substituting (15) and(16) into (14), we obtain 
31h? 8h4 
yy. a= 63 (3) (§2 7. 2 563) (63) (6M. 
5 DU; = 80; + 5) (6@U;) + 82 (6 (6OU;)) +, 


a3u, = 52U; —e 5(3) (62u,) - : Sh 3 (6 (6@u;)) —T. (18) 


Now substituting (17) and (18) into (2), we obtain semi-discretization of 
the given governing PDE and the obtained nonlinear system of ordinary 
differential equation (ODE) is 


dU; h? ht 
aU; (zy. (3)q7. 
a tus (si U; — —6QU; - 82 (08 U;) tr) 


v (se, - 4 ae aoe) (6205) - S52 ) (6 (6@u;)) - Ts) =0. 
This implies that 


dU; h? (3) ty h4 3) (2) ry ry 
= j dcx U; 4 Ocx (8 j) — Oca Uj 1 
a Us ( 6 7" 120 ) z oe) 


IJNAO, Vol. 12, No. 2, pp 277-300 


284 K. Aliyi Koroche and H. Muleta Chemeda 


31h? 53) ( 5(2) Bh* 03) (53) (5) (2) 
+uv (Seas (6&: U;) =F 315 Ge (65: U;)) | dcx U; ae ee 


5040 * 72576 * 
truncation errors with order of accuracy O(h°). Hence by truncating this 
local truncation error, our proposed scheme is 


dU; h? ht 
Be ( = dee Uj + bee (62u;) = 2205) 


1 
where T; = ( ges Or+ a 22) h°U;, is the jth terms of principal local 


dt 12 

+v (eae (52u;) + SEs (a (8@u;)) - ou; ) 
dU; _ah?U; (“2 QUj41 + 2Uj-1 vi-2) 
dt 6 2h3 

; vh4 (“ 4Uj42 + 5Uj41 — 5U;-1 + 4U;_2 “i= ) 

" 120 2h3 


apeeeas ia 
yf Ue = ier 
2h 


ve (es AU 542 + 5Uj41 — 5U;—-1 + 4U;~2 “i=t) 


20 2h3 
Gy ( Vite = 2541 + 2Uj-1 — Uj-2 
2h3 
128vh* / 1 
t 340 (= (Uj45 — 4Uj44 + 3Uj43 + 8Uj42 — 14U;41 


+14U;_1 — 4U;_2 — 3U;_3 + 4U;_4 Ui+s)). 


This implies that 


dU; U; 

dt a (Uj43 + 36Uj42 — 315U;41 + 315U;_1 — 36U;—2 — Uj—3) 
VU 

+ Grane (Uits — 4(Ui+4 — Uj—a) + 2607(Uj+3 — Uj—s)) 


+ (—30568(U; 42 — Uj—2) + 53326(U;+1 — Uj—1) — Uj—s), 


which implies that 


= a OE a8 RU NBO hes ie I ie cca) 
+ gggnn (249 — Uses) 
= sags (2(Uit4 — Uj—a) + 1803.5(Uj 49 — U;-4)) 
~ aggone (15284(Ui+2 - U;_2) + 26663(U;41 — Uj_-1)), 


(20) 


where j = 0,1,2,3,..., M—1. From (20), the nonlinear system of the ODE 
is given in the form of 
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dUy 3a 
Topp, (Obs — U-s) + 36(U2 — U_2) — 315(U; — 315U_1)) 
0.5u 
+ 3360n3 (Us — U-s) 
U 
33g43 (2(U4 — U-a) + 1303.5(Ug — U_s) — 15284(U2 — U_2) 


+ 26663(U; — U_1)) 


dU, _3aU 
“He = 1gop (0:1 U4 — U-2) + 36(Us — U1) — 315(U2 — 315U_o)) 
0.5u 
3360n3 (V6 ~ U-a) 
UV 
a360n3 (2(Us — U-a) + 1803.5(Us — U_2) — 15284(Us — U_1) 


+ 26663(U, — U_o)) 


dU 3aU ym 
gop eM a+4 — Um—s) + 36(Uae4s — Uar—1) 


— 315(Um42 — 315Um)) 
0.5u 


* 3360h3 (Um+6 — Um-a) 
1 


~ 39 gpq5 (2(Um+s — Uma) + 1303.5 + 4 — Uma) 


+ saqpq3 (15284(Um+s — U1) — 26663(Um+1 —Um)). (21) 


Now by introducing the vectors U(t) = [Ui(t), U2(t), Us(t),...,Unc(t)]’ in 


(21), we can rewrite the system of equation in the matrix form given by 


dU(t) 
— = AH U®) (22) 


subject to discretized initial and nondiscretized boundary conditions given 
by 
U(xj,0) =Upo(aj), OS 2; <b, 


U(0,t)=Up(t), U(b,0)=U,(t), O<t<T, (23) 


and where # is a nonlinear function of U which contains the entries h; give 
by H(U;) _ hj (Vi, U2, sees Uy-1). 
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2.2 Temporary discretization 


Consider the discretization of time domain [0,7] by using the equidistant 
mesh-size with time step At given as 


D, = Abr =tyht+ At} , (24) 


where At = T/N and n = 1(1)N, in which N is the Maximum number of 
grid points in temporary (time) direction. Then the resulting system of Odes 
in (22) can now be solved by using the classical fourth-order Runge-Kutta 
method. This Runge-Kutta techniques integrate the value of U from n to 
n-+1 by using the operator given by 


Uo = U(a;,t0); ky = AtAH(Up), 
Uy =Up+1/2Ky, Ky = AtAH (Ui), 
ih= y+, Ko= Aran), 
Us =Up+1/2Ks, Ka, = AtAH (U2), 
Unsi = Un + 1/6 (Ky + 2K2 4+ 2K3 4+ Ky), (25) 


with both discretized initial and boundary conditions given by 
Ue, 0) =Uye;), ON; <4, 


U(0, tn) = Uo(tn), U(b,0) = Uo(tn), O0<t, <T, 
where j = 1(1)M—landn=1(1)N -1. 


3 Stability, consistency, and convergence method 


3.1 Stability of scheme 


The stability of the proposed numerical method is investigated by using von 
Neumann stability analysis. To do this, we assume that the nonlinear term 
UU, of partial differential equation in (2) as linear by taking y = U;, where 
y is constant. Then without loss of generality, we obtain the linear system 
of ODEs. Assume that y = max,;(U;) in (20), and we can inquire about 
the eigenvalues of the M by M system of ODEs in (20). To obtain this 
eigenvalue, as [3, 12, 22, 25] takes, we assume that a trial solution is 


U(a,t) = p(t)e(a). (26) 


Furthermore to investigate the stability of the proposed method by von Neu- 
mann techniques, we assume that the trial function is defined by 
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(2) = eiKpx; =~ ee, (27) 
where 7 = /—1 and Kp = Pa, P = 1,2,3,...,M. Moreover, K is a Fourier 


number or amplification factor. Now substituting (26) and (27) into (20), we 
obtain 


dp(t)eIhKp 
dt 
= sre) (0.1(€i@XeG+8) — ethKp(G-3)) 4 g(ethkp G42) _ gee) 
7 sata (15(eh rot) _ eked) 
gee (et? Kp +5) = eth Kp(5—5)) 

3360h 

2 vat (2(e@oe+9 — eh Kp G-4)) 4 1303.5 (eth Kn G48) — eM>G-9))) 
3360h 

_ wee, (15284(ci™o+2) — cihkpG-2)) 4 26663(et* Ar +1) — eitkrG-1))) 


Dividing both sides by e4"*», this implies that 


dt 160h 
_ 8ayye(t) (15(cihK» — e-thKny) 4 0.5v (cdthKy _ @SihKy) 


dp(t) _ 3ayp(t) (0.1(c3"K — e-BIbKy) 4 4(e2thKy _ g-2thKp)) 


160h 3360h3 
vp(t) in ie: a 
ro 3360h3 (2(e4 hKp —e 4ihK, ) Fils 1303.5(e? hKp 6 3 hKp)) 
vp (t) 2ihK, —2ihK, thK. —ihK 
— aaqpns (15284(e™*Me — e BEM) + 26663(ePM” — er) 28) 


Since for any differentiable function y(t), the Eigen-value problem of y(t) is 


dy(t)/dt = Apy(t) 


for p= 1(1)M. Hence substituting this into (28), we obtain 


p(t) aetna (0.1(e8*Kp — e-3ihKp) 4 4(e24Ky _ e-2ihKp)) 
t F : 
_ aa (15(e**Kr _ gy) 
0.50 sink, —BihKp 
+ 336008 | _ ) 
t 
= a (2(etKr — e-4hKr) 4 1303,5(€8hKr — e~3ihK)) 
_ vet) 
3360h3 


(15284(e7""*» oz es) a 26663(e""*» _ ev) ; 


Therefore this implies that 
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t . ; 
dp _ 8ayo( ) (0.1(e3*X> — ¢ 8ihK) 4 4(e2ihKn _ e2thKy)) 


160h 
= a (15(e!*Ke — e~*hKr)) 4 oe (cdihKy _ @-SihKy) 
_ 3360K3 Ce _ ona a 1303.5(€34Kp = ee) 
_ 3360K3 (15284(e?*"Ke — e~2hKe) 4 96663(eiKe — e-hKy)) (29) 


From (29), using the Euler expiration into trigonometric form, we obtain 


on =a (0.5 sin(5hKp) — 2sin(4hK,) + (6.3yh? + 1303.5) sin(3hKp)) 
252yh? — 15284) sin(2hKp) + (945 yh? — 26663) sin(hKp)) . 


(30) 


‘eae 
1680h3 


Hence from (30), we obtain the required Eigen-value. Real part of all Eigen- 
values is zero. This shows that strictly Real(\,) < 0 and so that, we have 
|Ap| < 1 for all values of variable ”p”. Thus the required condition is satisfied. 
Therefore the obtained system of equation in (20) is stable. 


Theorem 1. The difference mathematical statement given in (20) is sta- 
bilized equation if and only if all eigenvalues of coefficient matrix of ob- 
tained system of difference equation are simultaneously satisfy condition 
Real(Aq) < 0. 


Proof. The proof of this theorem is briefly given in [22]. 


3.2 Consistency of the scheme 


Round-off errors and truncation errors occur when mathematical equations 
are solved numerically. Rounding errors originate from the fact that com- 
puters can only represent numbers (approximate solution) by using a fixed 
and limited number of significant figures. Therefore, such numbers cannot 
be represented exactly in different computer memories. The error introduced 
by this limitation is called a round-off error. Truncation errors in numerical 
investigation arise when approximations are used to estimate by using the 
finite term for infinite quantity by truncating certain terms from the bound- 
less series term. Note that in numerical computation, if we work to minimize 
round-off error, then the truncation error is increased. Therefore it is diffi- 
cult to avoid these errors in numerical computation. So that, in our work, 
we must find the balancing (equilibrium) points and minimize both of them 
simultaneously. 

Generally, in all proposed numerical difference schemes, the accuracy of 
solution for a given mathematical problem obtained by those methods de- 
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pends on how small we make the step size, h, and time step At. Therefore, if 
the local truncation error produced in a numerical scheme is near to zero for 
which both h and At simultaneously approach to zero, then this numerical 
method is called consistent. 

Hence we have in-minding this, from our present proposed scheme in (20), 
that the principal part of local truncation error (the derived local shortening 
error) is 


lim T; = lim ( 


j-co jroo 


U; 51 
o* + —__ 8° | A°U; = 0. 
5040°" | 12096 ") é 

Therefore from the above equation and the proposed theoretical-bounds of 
truncation error, we have max;||Zj|| < O (h°). Hence this implies that 
max; ||Z;|| > 0 as h + 0. So that, our scheme is consistent with the order 
of accuracy O (h®). 


3.3 Convergence analysis 


Stability plus consistency of the proposed scheme leads to the convergence of 
the method. Hence our present system is both consistent and stable. So that, 
the outline present method is convergent with six-order of convergence in the 
spatial direction. To test the performance of this convergent and accuracy 
of the proposed method, we use the maximum point-wise absolute error, Lz 
and L,, norms. Theses norms are calculated by 

Loo = max; |U(2;,tn) — Ujn| and Ly = /1/M Se GU al 
where M is the maximum number of grid points, U(t;,t,) is exact and Uj, 
is an approximation solution of the KdV-Burger equation in (2) at the grid 
point (x;,t,). If the exact solution of the model problem does not exist, then 
we take the approximate solution U2; 2, © U(x;,tn) obtained at (x;,tn) by 
taking double number of grid-point as an exact result, and we can compare 
the error between them at specified grid points. 


4 Results of numerical experiments 


To test the validity of proposed method, we have considered the following 
model problem. 


Example 1. Consider that one-dimensional KdV equation given by 
U,+a(U*) + BU rac = 0 
with periodic boundary conditions and the following initial condition: 


(A) Uo(a) =3 a sech(G(x — c)) +3 b sech?(P(x — d)) 
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a=0.3, b=0.1, c=0.5, 0< 2 <2,0<t<4, 
B = 5.84 x 1074,a = 0.5,G = av/a/8, P = av/b/f, 
(B) Uo(x) =3 a sech(G(a — c)), 
a=0.38, c=0.5,0<¢4<1,0<t<4, 
B=10-*, G=1/2V/a/8, 


2 (az — 1) 
CO) Ile) =— h? 
(C) Vola) = 5 a sech*( =), 
a=0.5,0<27<3,0<t<4, 


B=10%, 


Example 2. Consider the one-dimensional KdV equation considered in [37] 
given by 
U, + 6UUz + Urexr = 0. 


With analytical solution obtained in [37] given by 


U(z,t) = 5 sech”( (x — ct — D)), 


2 


where D is a constant of integration. 
Example 3. Consider the KdV equation considered in [30] given by 
U, + (a+ BU)UU, + Wer — Uzre = 9, 


where initial condition of this equation for ~ = 0 and 6 = —1 is given by 


U(«,0) = (1 + tanh (<3) 


In [80], by letting 6 = 0, a = 2, y = —5, and 6 = —3, its analytical solution 
is given by 


U(az,t) = : (sexn*(5) + 2 tanh(5) +2) ; 


where 6 = $2 + 3t. 


Example 4. Consider that fractional KdV equation considered in [1] given 
by 
O°U OU &U 
+ U = F(z,t), (x,t) € (0,1). 
ote Ox 


Analytical solution obtained in [1] is given by 


U(a,t) = te”, 
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by assuming a = 0.9. 
Table 1: Point-wise absolute error and root mean square error, for Example 1 (A) if 


M =61 
provided grid points estimated error at provided grid point 


x t Leo Ly 
0.66667 3.9333 4.7839. — 03 1.1526E — 04 
1.6667 3.9833 4.4147. — 05 3.11654 — 06 
1.9667 3.9983 3.5070. — 05 2.300E — 06 

2 4 2.4568.E — 06 2.300E — 06 


propagation at t=O propagation at t=1 


Behaver of solution 


Figure 1: Graph of the physical behavior of the solution when 6 = 4.25 x 10~-°4 and 
M = 61 with high-frequency oscillations for Example 1 (A). 


5 Discussions 


In this paper, sixth-order compact finite difference method is presented for 
solving a one-dimensional KdV-Burger equation. At first stage, we discretized 
the solution interval for concerning spatial variable. Next by applying the 
Taylor series expansion, we discretization partial derivative of the model prob- 
lem in (2) involving spatial variable, we obtain sixth-order finite difference 
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Table 2: Point-wise absolute error and root mean square error, for Example 1 (B) if 
M=61 
provided grid points estimated error at provided grid point 


x t Leo Lo 
0.333 3.967 4.7 — 04 1.3838 — 05 
0.8333 3.9917 3.600 — 05 32.2569E — 06 
0.9833 3.9992 3.532E — 06 1.3216 — 07 
1 4 2.473E — 07 3.2814 — 08 


Table 3: Point-wise absolute error and root mean square error, for Example 1 (C) if 
M=61 
provided grid points estimated error at provided grid point 


x t Lixo Lo 
1.333 3.8667 7.2E — 04 3.6786. — 05 
3.333 3.9677 9.7 — 05 4.0591F — 06 
3.9333 3.9967 6.25E — 07 5.1082E — 07 

4 4 3.348 E — 07 3.126E — 08 


Table 4: Point-wise absolute error and root mean square error obtained by the present 
method, for Example 2 using different step sizes h and time step 
provided step size and timestep Estimated Error at provided grid point 


h At Leo Lz 

0.1 0.5 8.128E — 04 2.1526E — 05 
0.05 0.5 5.9148E — 05 7.6374E — 06 
0.02 0.1 4.1139E — 05 5.2931 E — 06 
0.01 0.1 3.8425E — 06 9.2631E — 07 


Table 5: Point-wise absolute error and root mean square error obtained by the present 
method, for Example 3 with different step sizes h and time step 
provided dance step size and time step Estimated Error at provided grid point 


h At bse L2 

0.1 0.04 1.3298E — 06 6.3691E — 07 
0.1 0.03 5.1653. — 07 4.6334E — 07 
0.1 0.02 4.623E — 07 4.211E — 07 
0.1 0.01 1.975E — 07 8.091F — 08 
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propagation at t=O propagation at t=1 


propagation at t=2 Behaver of solution 


Figure 2: Graph of the physical behavior of the solution when 6 = 4 x 107-4 and 
M = 61 with high-frequency oscillations for Example 1 (B). 


propagation at t=0 propagation at t=1 

1 1 

08 08 
06 06 
$04 $04 
0.2 0.2 

0 0 

0 1 3 4 0 1 3 4 
x x 


propagation at t=2 


Figure 3: Graph of the physical behavior of the solution when 6 = 10-°4 and M = 61 
with high-frequency oscillations for Example 1 (C). 


IJNAO, Vol. 12, No. 2, pp 277-300 


294 K. Aliyi Koroche and H. Muleta Chemeda 


a) Surface of exact solution 


time 0 -10 x axis 


b) Above-ground of approximation result 


a 
05 ~ oe ) 


time 0 -10 x axis 


Figure 4: Graph of the physical behavior of the surface of the present numerical solution 
with mesh sizes h = 0.01 and k = 0.01 for Example 2 which we compare it with the 
surface of exact solution in [37] 


scheme. Then, by rearranging these semi-discretized schemes and combin- 
ing them into the remaining part of partial derivative in the model problem 
concerning temporal variable without linearizing it. We obtain a system of 
nonlinear ODEs. Then, we solve the obtained system by using the fourth- 
order Runge-Kutta method. To demonstrate the competence of playacting 
(to identify the applicability and validity of the present proposed scheme), 
four model examples are solved by taking different values for step sizes h 
and dispersion coefficients v, a, and @ in the different solution domains. To 
further verify this validity of the proposed method, the norm of errors be- 
tween the numerical results for model problem in Example 1 subjected to 
three different initial conditions associated with numerical results obtained 
by dabbling mesh size (2M,2N) taking as the exact solution, is summarized 
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a) graphs 
| 


0.9 
08 
07 


| 0 
x 


c} surfece of Numerical solution for example three 


of numerical versus exact solution of Example three 


—— u-approximat for t=0.9 
—- U-exact for t=0.9 


| 


1. 
b) surfece of Exact solution for example three 


d) Physical bahavertof emror Belween exact Versus Numerical 
solution for example three 


u(x,t) 


6 


x10 


295 


Figure 5: Graph of the physical behavior of the surface of the present numerical solution 
with mesh sizes h = 0.05 and k = 0.01 for Example 3, which we compare it with the 
solution in [30] 


Table 6: Comparison of point-wise absolute error and root mean square error, for 
Example 4 using different step sizes h and time step At 


provided step size and time step 
result by present method 


h At 
1/10 1/20 
1/30 1/20 
1/50 1/20 

Result by Method in [1] 

1/10 1/20 
1/30 1/20 
1/50 1/20 


Estimated Error at 


L 


4.5349. — 05 
5.1653 — 06 
3.5283. — 06 
1.5678E — 04 
7.8531E — 05 
1.6579E — 05 


Ly 
8.3691E — 06 
5.6334E — 06 
9.4184E — 07 
6.1032E — 04 
3.9868E — 05 
7.1991E — 06 


provided grid point 


in Tables 1-3 and Figures 1-3 for different disparate parameters and its so- 
lution domain. 

As it can be seen from these tables, the point-wise maximum absolute 
errors (L..) and root mean square error (£2) are rapidly decreases as the 
value of x and t increase. This indicates that the vibration of the solution 
line has to become overlapping and that the dispersion of the travelling wave 
is equal. From those figures, also we can understand that, for different values 
of 6 and initial conditions, we obtain differences in diametrical travelling 
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a) graphs of numerical versus exact solution of Exaple four 
T T T T 1 


—}— w-approximat for t=1 
wescene —+ U-exact for t=1 ee. Semen 


ucx.t) 


Maximum Absolute error 


grid poit of in spatial direction 


Figure 6: Graph of the physical behavior and surface of the present numerical solution 
with mesh sizes h = 0.01 and At = 0.1 for Example 4, which we compare it with the 
surface of exact solution in [1] 


wave, and at a certain point in the assorted domain, the vibration line of the 
travelling wave is overlapping. Hence, in this case, a numerical approximation 
of Example 1 for different values of 8, initial conditions, and the solution 
interval is very well, and it gives different behaviors of the solution. From 
Table 4 and Figure 4, we conclude that the present numerical resolution for 
Example 2 is accurate than the quantitative result in [37]. Thus this not only 
shows the solution given by the graphical form, but also as we saw from Table 
4, if mesh-size decreases, both roots mean square error and point-wise error 
of the present numerical result are rapidly decreased. So that, this shows 
that the accuracy of the present acting method is increased and proves the 
preexisting accuracy. Hence the approximate solution by the present method 
for Example 2 is very well. Also, as we saw from Table 5 and Figure 5 that 
the present numerical solution of Example 3 is accurate than the numeric 
resolution in [30]. Because of, as it can be seen from the obtained error 
norm of numerical results specified in Table 5, both roots mean square error 
and maximum point-wise error are decreased. This shows the correcting and 
proving the preexisting results in [30]. Hence our result is more accurate 
than the event that exists in the literature. In addition to this, the present 
techniques are also used to solve a fractional order KdV-Burger equation. As 
it can be seen from the comparison of accuracy through the error norm for 
numerical results given in Table 6, the quality of the present method is more 
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superstar than the truth of solution in [1] for a = 0.9. Because of both Loo 
and Lz» rapidly decrease as h decreases to zero with a fixed time step At. 
Therefore, depending on all accuracy of the present numerical results, the 
proposed method validates and improves the accuracy of preexisting meth- 
ods described in the literature. Also in this paper, both the theoretical and 
numerical error bounds have been established. The obtained quantitative 
results listed in the tables and graphs are also confirmed our established com- 
putational rate of convergence and theoretical estimates of the error bound. 


6 Conclusion 


In this study, a new approach, which was called sixth-order compact finite 
difference method, was presented to solve a one-dimensional KdV-Burger 
equation. The comparison of results obtained by the present method acting 
was more convenient, reliable, and effective than some methods listed in the 
literature. An error analysis based on the von Neumann stability investigat- 
ing was also developed for the present study. So that, as it can be seen from 
both graphs and tables, the accuracy of the present method has improved 
the inaccuracy of preexisting techniques in literature by decreasing the value 
of h and At. As a summary, the accuracy of numerical results given in terms 
of tables and graphs indicated that as the values of h and At decrease, the 
errors drop-off more rapidly. Generally, this technique is a reliable acting 
method and ables solving the one-dimensional KdV-Burger equation. Based 
on the findings, this method is well approximate and gives better accuracy 
of the numerical solution with both decreasing or some fixed step sizes h and 
time step At. Therefore this method is a suitable technique to approximate 
the exact solution of PDE very well. 
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